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Abstract 

We use a hybrid N-body program to study the evo¬ 
lution of massive black hole binaries in the centers of 
galaxies, mainly to understand the factors affecting 
the binary eccentricity, the response of the galaxy to 
the binary merger, and the effect of loss-cone deple¬ 
tion on the merger time. The scattering experiments 
from paper I showed that the merger time is not sen¬ 
sitive to the eccentricity growth unless a binary forms 
with at least a moderate eccentricity. We find here 
that the eccentricity can become large under some 
conditions if a binary forms in a galaxy with a flat 
core or with a radial bias in its velocity distribution, 
especially if the dynamical friction is enhanced by res¬ 
onances as suggested by Rauch and Tremaine (1996). 
But the necessary conditions all seem unlikely, and our 
prediction from paper I remains unchanged: in most 
cases the eccentricity will start and remain small. As 
a binary hardens it ejects stars from the center of a 
galaxy, which may explain why large elliptical galax¬ 
ies have weaker density cusps than smaller galaxies. If 
so, the central velocity distributions in those galaxies 
should have strong tangential anisotropies. The wan¬ 
dering of a binary from the center of a galaxy miti¬ 
gates the problems associated with loss-cone depletion 
and helps the binary merge. 

1 Introduction 

The formation of massive black hole (bh) binaries fol¬ 
lows from two widely held assumptions: that galax¬ 
ies merge and that many galaxies contain central bhs. 
The evolution of these binaries is relevant to a number 
of problems in extragalactic astronomy. They have 
been proposed to explain the bends in jets from ac¬ 
tive galaxies (Begelman et ah, 1980; Roos, 1988; Roos 
et ah, 1993; Gaskell, 1996), and the difference between 
radio-loud and radio-quiet active galaxies (Wilson and 


Colbert, 1995). Multiple mergers of bh binaries can 
lead to the buildup of massive central bhs (Hut and 
Rees, 1992), although if the binaries do not merge 
quickly the bhs can be ejected through three-body in¬ 
teractions (Xu and Ostriker, 1994; Valtonen, 1996), a 
process that has been proposed to explain the double 
radio lobes of active galaxies (Saslaw et ah, 1974; Val¬ 
tonen et al. 1994). The double nucleus of the nearby 
galaxy M31 (Lauer et al., 1993) could be a bh bi¬ 
nary in the making, although Tremaine (1995) has 
proposed a plausible model that requires only one bh. 
If binary mergers are frequent they will be the most 
interesting source of gravitational radiation for the 
planned space-based detector LISA (Haehnelt, 1994; 
Bender et ah, 1995). Lastly, the evolution of bh bi¬ 
naries may explain why the density profiles of ellipti¬ 
cal galaxies fall into two classes, with small galaxies 
having strong density cusps and large galaxies hav¬ 
ing weaker cusps (Ebisuzaki et al., 1991; Faber et al., 
1996). 

The relevance of BH binaries to all of these problems 
depends on how long it takes for dynamical friction 
(resulting from star-BH interactions) and gravitational 
radiation to cause the binaries to merge. The merger 
time is difficult to compute because the hardening of 
the binary changes the structure of the galaxy, which 
in turn changes the hardening rate. The back-of- 
the-envelope calculations of Begelman et al. (1980) 
have now been supplemented with numerical results 
from three-body scattering experiments (Roos, 1981; 
Mikkola and Valtonen, 1992) and full N-body exper¬ 
iments (e.g., Ebisuzaki et al., 1991; Makino, 1997), 
yet several big uncertainties remain. The eccentricity 
evolution has been subject to much debate, the most 
recent contribution being that of Rauch and Tremaine 
(1996) who argue that resonant interactions between 
the BHS and the stars can cause a rapid growth in 
the eccentricity if the velocity distribution has a ra¬ 
dial anisotropy. This would cause gravitational radi- 
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ation to become important early in the evolution and 
would help the binaries merge. Another uncertainty 
has been loss-cone depletion, the depletion of the low- 
angular-momentum stars that can interact with the 
binary, which lowers the hardening rate and length¬ 
ens the merger time. Because of the uncertainties one 
can find papers in the recent literature making widely 
different predictions for merger times in typical galax¬ 
ies, ranging from <10® yr (e.g., Makino and Ebisuzaki, 
1994) to well over a Hubble time (e.g., Rajagopal and 
Romani, 1995; Valtonen, 1996). 

In paper I (Quinlan, 1996) we presented results 
from a new set of three-body scattering experiments 
as a first step towards resolving the uncertainties. We 
computed the hardening rate, the eccentricity growth 
rate, and the rate at which a binary ejects stars, and 
derived formulas to fit the dependence of the rates on 
the hardness and mass ratio of the binary. One of the 
main results was that the eccentricity of a hard binary 
grows only slowly as the binary hardens, and that the 
growth does not change the merger time by much un¬ 
less the binary starts with at least a moderate eccen¬ 
tricity, say e > 0.5 (this had been demonstrated for 
equal-mass binaries by Mikkola and Valtonen, 1992). 
But in deriving these results we made some simpli¬ 
fying assumptions that may not be applicable to real 
galaxies: we assumed a homogeneous, isotropic galaxy 
core with a Maxwellian velocity distribution, and we 
ignored the change in the core caused by the bi¬ 
nary. We thus could not study the response of the 
galaxy and the depletion of loss-cone orbits in a self- 
consistent manner, nor could we study the dependence 
on the initial density profile or on anisotropies in the 
velocity distribution. 

To improve upon the scattering experiments we 
have developed a hybrid N-body program that is well 
suited to studying the evolution of bh binaries in mod¬ 
els for elliptical galaxies, including models with cusps 
and anisotropies. In this paper we describe the de¬ 
sign and performance of the program and then use it 
to study the factors affecting the eccentricity, the re¬ 
sponse of the galaxy to the binary hardening, and the 
effect of loss-cone depletion on the hardening rate. We 
pay special attention to the convergence of the results 
with the number of particles and to some differences 
between the N-body results and the predictions from 
paper I. 

We ignore some complications that may be impor¬ 
tant for real galaxies. We consider only spherical or 
nearly spherical galaxies, and ignore disks or triaxial 
components of the potential. A disk will alter the fric¬ 
tional force on the bhs, possibly dragging them into 
the disk plane, and a triaxial component will help the 
merger by increasing the number of stars that can 
interact with the binary. We also ignore gas, which 


can help the merger by accreting onto the bhs or by 
applying a frictional force to the bhs. We wish to 
understand the simple problem with purely stellar- 
dynamical processes first before we add complications. 
The neglected complications should not matter for bi¬ 
naries in gas-poor, nearly spherical galaxies like M87. 

2 Computational methods 

2.1 The N-body program 

N-body programs are usually designed to solve prob¬ 
lems of one of two types: collisional or collisionless. In 
collisional problems the evolution depends in an essen¬ 
tial way on close encounters between the stars. These 
problems require accurate integrations with the exact 
1 /r^ force law, and are usually solved with high-order 
integrators. In collisionless problems the evolution de¬ 
pends not on close encounters between the stars but 
on the collective response of the system to changes in 
the potential. These problems require a large num¬ 
ber of particles to suppress spurious relaxation, and 
are usually solved with a low-order integrator (e.g., 
leapfrog) and a fast, approximate method for com¬ 
puting the forces. The BH-binary problem does not fit 
neatly into either of these types. The binary evolves 
because of close encounters between the bhs and the 
stars, which must be integrated accurately, but the 
galaxy evolves because of the collective response of 
the stars to changes in the potential, and the number 
of stars must be large, both to suppress relaxation and 
to satisfy the requirement that the stars be much less 
massive than the bhs. 

It is not surprising then that both types of N-body 
programs have been used for the problem. Governato 
et al. (1994) used a tree program with the leapfrog 
integrator and a large softening length to study the 
merger of King models containing central bhs, with 
16,000 particles per galaxy. Their program was good 
for the galaxy merger and the early stages of the bi¬ 
nary evolution, but not for the late stages, which they 
did not study. Mikkola and Valtonen (1992) used 
Aarseth’s NBODYl program to study the evolution 
of a BH binary at the center of a Plummer model with 
10,000 particles. NBODYl uses a fourth-order inte¬ 
grator with individual stepsizes, which makes it better 
than a tree program for the late stages of the evolu¬ 
tion, but it computes the forces by direct summation, 
which is slow. Mikkola and Valtonen assumed a fixed 
Plummer potential when computing the forces be¬ 
tween the stars, which made the integrations go faster 
but at the price of ignoring the self-consistent response 
of the galaxy. Self-consistent integrations of galaxy 
mergers with central bhs using over 10® particles per 
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galaxy were done with a direct-summation program 
by Makino and Ebisuzaki (1996) and Makino (1997), 
but on a GRAPE-4 supercomputer, whose speed of 
100-150 Gflops is beyond the reach of general-purpose 
computers. 

Our hybrid N-body program works well for both the 
collisional and collisionless aspects of the BH-binary 
problem. We call it SCFBDY because it combines 
parts from the SCF program of Hernquist and Os- 
triker (1992) and the NBODY programs of Aarseth 
(1994). We describe here the main features of the pro¬ 
gram and their advantages for the BH-binary problem; 
a more technical description is given in the appendix. 

The program moves the bhs and stars according to 
different equations of motion (we use lower-case letters 
for the stars, upper-case for the bhs, and set G = 1): 

( 1 ) 


rk = Yl 
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The BHS interact with the other bhs and with the 
stars through the 1 /r^ force (a softening parameter 
is added for the BH-star interactions). The stars in¬ 
teract with the BHS through the 1/r^ force, but with 
the other stars through an expansion of the poten¬ 
tial with coefficients Anim that are updated with time 
self-consistently. By using the ^nim basis functions 
of Hernquist and Ostriker (1992) we can fit the po¬ 
tentials of most elliptical-galaxy models adequately 
with a small number of coefficients—a dozen is often 
enough for spherical galaxies. Phinney and Villumsen 
(Phinney 1994) split the equations for the BH-binary 
problem in a manner similar to our equations (|^ and 
(^, but used a different expansion for the potential 
and used the leapfrog integrator to move the parti¬ 
cles. We move the particles with individual stepsizes 
(sometimes differing by as much as a factor of 10 ®) us¬ 
ing one of Aarseth’s (1994) fourth-order integrators— 
the NBODYl integrator for the stars, and NBODY2 
or NBODY 6 for the bhs. Our program thus combines 
the integration accuracy of a direct-summation pro¬ 
gram with the speed of a SCF type of program. 

The program does have some limitations. It does 
not vectorize or parallelize as well as the pure SCF 
program. It modifies two-body interactions between 
the stars by replacing the potential by a smooth ap¬ 
proximation, and hence cannot be used to study the 
diffusion of stars into the loss cone by two-body re¬ 
laxation. (Programs that do not replace the potential 
by a smooth approximation commit the opposite er¬ 
ror for the BH-binary problem: they give a diffusion 


rate larger than would occur in a real galaxy, because 
they cannot use as many stars as there are in a real 
galaxy.) And the program requires a central point 
about which the potential can be expanded, at least 
in the simple version used here with only one expan¬ 
sion center, which means that it cannot be used to 
study the galaxy mergers that lead to galaxies with 
more than one central bh. We ignore that stage of 
the evolution. 


2.2 Initial conditions and nnits of mea¬ 
surement 


We used three simple, isotropic galaxy models for our 
experiments in this paper, as well as some anisotropic 
models that will be described later. The first of the 
simple models is the Plummer model, with a density 
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The Plummer model is easy to set up and interpret, 
but is a poor model for real galaxies: they are much 
more concentrated, many with densities that continue 
rising to the innermost observable radius (Crane et ah, 
1993; Ferrarese et ah, 1994; Forbes et ah, 1995; Lauer 
et ah, 1995). We therefore used two models with 
density cusps from the family considered by Dehnen 
(1993) and Tremaine et al. (1994): 


3-7 Md 
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We used the models with 7 = 1 and 7 = 2 , which 
are the models of Hernquist (1990) and Jaffe (1983). 
The “7 models” fit density profiles of elliptical galax¬ 
ies reasonably well if the length scale d is chosen to 
be close to the effective radius. The inner parts of 
elliptical galaxies can be fit better by Zhao’s (1996) 
formula, p{r) ~ with (3 typically 

in the range 2-3 (the length scale d, or “break radius”, 
is much smaller than the effective radius in these fits), 
which is similar to the fitting formula used by Byun et 
al. (1996). The difference between Zhao’s models and 
the 7 models is unimportant for our work because 
we are mainly interested in the region r < d where 
the models are the same. We truncate the models at 
r = 300d. 

To increase the statistical resolution near the center, 
and to help satisfy the requirement that the stars be 
much less massive than the bhs, we give the stars a 
mass spectrum, with the masses decreasing towards 
the center as 

m~r^, (5) 

with Vp the initial pericenter of the star and A an ex¬ 
ponent, usually in the range 0.5-1.0 (Sigurdsson et 
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al., 1995). In a Jaffe model with A = 0.75, for exam¬ 
ple, about 40% of the stars have r/d < 0.1 and 20% 
have r/d < 0 . 01 , whereas the corresponding percent¬ 
ages for an equal-mass model are 10% and 1%. The 
mass spectrum does not lead to mass segregation, be¬ 
cause the potential is smooth with our program and 
the segregation time is long. We impose minimum 
and maximum masses, typically differing by a factor 
of 1000 , so that the masses do not become arbitrarily 
small or large at small and large radii. 

The BHs are treated as point masses obeying Newto¬ 
nian gravity. The detected massive bhs in real galax¬ 
ies have a mean mass of 0.002-0.003 (in units of the 
bulge mass), although there is a scatter of at least one 
order of magnitude about this mean (Kormendy and 
Richstone, 1995). It is difficult to use small bhs in 
the experiments because the dynamical-friction time 
for those is long and because we need the stars to 
be much less massive than the bhs. We therefore 
use BHS somewhat larger than occur in real galaxies, 
'mhh/M ~ 0 . 01 , and vary their masses to determine 
how the results scale with the mass. The scaling is 
not important for the eccentricity results, which are 
not sensitive to the masses (except for the results on 
resonant friction), but it is important for the response 
of a galaxy to the binary evolution. Our masses are 
not that unrealistic given that our galaxy models have 
less mass at large radii than a real galaxy would have 
(if we interpret r = d as the break radius). 

The initial condition for our bhs are of two types. 
In the first we start the bhs symmetrically about the 
origin in one of the galaxy models described above, at 
r ~ d. We assume that the bhs arrived there after a 
galaxy merger, or after being dragged in from the halo 
of the galaxy. After a merger the bhs would proba¬ 
bly be surrounded by clusters of bound stars, but we 
ignore those (our potential expansion has trouble re¬ 
solving small clumps of stars away from the center). 
In the second we start a galaxy with one large bh at 
the center, surrounded by a cluster of bound stars, 
and allow a smaller bh with no cluster of bound stars 
to sink in from r cs d. The second type of initial con¬ 
dition is more realistic than the first, but the first is 
adequate for answering some questions and is easier to 
use than the second, because in the second the central 
BH and the stars surrounding it require small stepsizes 
right from the start of the integration. 

We present the N-body results in a system of units 
in which G = M = 1, with M the mass of the galaxy 
without the bhs. We usually choose the length scale so 
that E = —1/4, with E the energy of the initial galaxy 
without the the bhs. The length scales d for the Plum¬ 
mer and 7 models are then Stt/IG and 1/(5 — 27 ), and 
the half-mass radius and the density and dynamical 
time at that radius are all of order unity. We specify 


the softening for the BH-star interactions by either the 
softening length e or the softening velocity 

Ve = \/Gmh\i/e. ( 6 ) 

As in paper I, we define the binary orbital velocity 
Vbin by 

Vbin = x/CMia/a, (7) 

with Mi 2 = mi -|- m 2 . 

3 A sample integration 

We shall describe one N-body experiment in detail to 
explain our assumptions, integration procedure, and 
method of data analysis, and to identify some of the 
questions to be studied later. We also use this ex¬ 
periment to test the sensitivity of the results to the 
integration parameters and the number of particles. 

The problem that we consider is the hardening of 
an equal-mass binary at the center of a Jaffe model. 
The galaxy was represented by 10® particles with 
a mass exponent A = 0.75. The bhs had masses 
mi = m 2 = 0.01, and were started on nearly cir¬ 
cular orbits at (r, u) and (—r, —u), with r = 0.5. 
We integrated the orbits for 80 time units with an 
accuracy parameter rj = 0.005, a softening velocity 
Ve = 16, and a spherical potential expansion (with 
Umax = 18, l-max = 0) whose coefHcients were updated 
every At = 0.005; the integration took five days of 
cpu time on a DEC 250-4/266 AlphaStation. The 
program regularized the orbits at time t = 17.85. The 
total energy was conserved to about one part in 100. 

The integration results are shown in Figure |^. The 
semimajor axis a and eccentricity e were computed 
from the bh coordinates and velocities by the usual 
formulas for a Kepler orbit, ignoring the interaction 
between the bhs and the stars. The eccentricity re¬ 
mained small (e;< 0 . 1 ) and is not shown in the figure. 
The semimajor axis is shown in panel (b). The pas¬ 
sage of 1 /a through zero marks the time when the bhs 
become bound, t = 6.15, at which the separation be¬ 
tween the BHS is r = 0.063. A bh binary becomes hard 
once its semimajor axis shrinks to Oh = G'm 2 / 4 cr^ 
(Quinlan, 1996), with a the one-dimensional disper¬ 
sion near the center but outside the Keplerian rise 
around the bhs. For the Jaffe model, l/oh = 200 
(ct^ = 1/2). The hardening slows down noticeably af¬ 
ter 1/a reaches about 500, because the binary has by 
then ejected many stars and has reduced the central 
density. But the hardening does not stop, because 
the binary does not remain at the center: it wanders 
around and is thus able to interact with new stars. 
This is apparent from the plot of the center-of-mass 
radius in panel (a). 
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Figure 1: Binary evolution in a Jaffe model. Two equal-mass bhs (mi = m 2 = 0.01) were started on nearly circular 
orbits at (r, v) and (— r, —v), with r = 0.5. The panels show (a) the separation between the bhs (the lower line towards 
the end) and between their center of mass and the origin, (b) the semimajor axis, and (c) the radii containing fixed 
percentages of the galaxy mass: 0.008 (bottom line), 0.04, 0.2, 1, 5, 20, 50, 80, 95, 99, 99.8, and 99.96 (top line). The 
dashed line in panel (b) is from an integration in which the center of mass of the binary was fixed at the origin. 


To check the importance of wandering we repeated 
the integration with a constraint force applied to keep 
the center of mass of the binary fixed at the origin. 
We did this by moving the bhs with the same step- 
sizes and by replacing the forces T) and F 2 on the 
BHS at each step by Fi — F^m and F 2 — Fcm, with 
^cm = i'miFi -I- m 2 F 2 )/Mi 2 . With this constraint 
force the hardening stopped after 1/a reached about 
650—see the dashed line in panel (b). A series of inte¬ 
grations like this with different binary masses showed 
that the hardening stops when the mass in stars that 
can approach within a of the center drops to a small 
fraction of the binary mass. For the Jaffe model the 
final, “loss-cone” separation aic therefore varies lin¬ 
early with the binary mass, ajc ~ M 12 . This differs 
from the prediction of Begelman et al. (1980), but is 
similar to the predictions made by Roos (1981). When 
the binary is free to wander the hardening does not 
stop; we shall say more about this later. 

We stopped the integration at t = 80 because it was 


becoming slow, owing to the small stepsizes of the bhs. 
In a real galaxy gravitational radiation would at some 
point cause the bhs to merge. The time at which 
that would happen in our experiment depends on the 
physical mass and length scales, which we have not 
specified. Our neglect of gravitational radiation is jus¬ 
tified during most of the evolution because the impor¬ 
tance of radiation turns on suddenly—the timescales 
for hardening by three-body scattering and by radia¬ 
tion vary as ~ 1/a and ~ a^. Our goal is to integrate 
the Newtonian equations far enough that we can pre¬ 
dict how the evolution would continue until the point 
at which radiation becomes important. For the appli¬ 
cations that we are interested in— bhs of mass 10 ®- 
IO^Mq in a large galaxy like M87, or bhs of mass 
10 ®-10 ® Mq in a smaller, denser galaxy like M32— 
radiation becomes important after the semimajor axis 
shrinks by a factor of 10-50 beyond a = ah (Quinlan, 
1996). In our integration the semimajor axis shrank 
by a factor of 12 beyond a = ah, which might have 
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been enough. 

Panel (c) shows how the galaxy expands as the stars 
absorb energy from the bhs. This starts when the bhs 
first sink into the center, before they become bound, 
but most of it occurs after the binary forms and begins 
hardening. The expansion appears to slow down near 
the end, but that is mainly because the figure uses a 
linear scale for 1/a and a logarithmic scale for r. Near 
the end the outer parts of the galaxy begin expanding 
too, because of the stars that the binary has ejected. 

To test the dependence of the results on the inte¬ 
gration parameters we repeated the integration with 
different choices of 77, At, nmax, and tmax- None 
of the changes had a significant effect on the binary 
evolution unless a parameter was changed to an unrea¬ 
sonable value, such as a softening velocity lower than 
the orbital velocity of the binary at the end of the in¬ 
tegration (T4in = 6.8), or a value of rj or At so large 
that energy conservation was grossly violated. Even 
adding non-spherical terms to the potential by chang¬ 
ing tinax from 0 to 2 hardly affected the binary, because 
the galaxy remained nearly spherical throughout the 
integration. The parameters that we used thus seem 
adequate for our purposes. 

To test the dependence on the number of parti¬ 
cles we compared results from integrations using 6250, 
25000, and 10® particles, using both equal-mass par¬ 
ticles (A = 0) and particles with a mass spectrum 
(A = 0.75). For each (A, A) combination we per¬ 
formed three integrations, with the orbits of the BHs 
started in the XY, YZ, and ZX planes; the difference 
between them is a measure of the noise level. The 
results are summarized in Figure 

The eccentricity is noisy when N is low, especially 
with the models with equal-mass particles, because 
the eccentricity of a hard binary gets perturbed only 
by a small fraction of the stars—those that have close 
encounters with the binary. In the models with a mass 
spectrum that fraction is larger and the eccentricity 
is less noisy. The convergence of the eccentricities 
towards e ~ 0.05 as N rises in panels (d), (e), and (f) 
suggests that it is the correct answer. 

The semimajor axis is not as noisy as the eccen¬ 
tricity. Even with only 6250 particles the integrations 
starting from the three orbital planes give nearly the 
same hardening rate. But there is a systematic de¬ 
crease in the rate as N rises from 6250 to 10®. A 
similar decrease was observed by Makino (1997), who 
proposed two explanations: two-body relaxation and 
the wandering of the binary from the center of the 
galaxy. Since our program suppresses two-body re¬ 
laxation, the decrease most likely results from inter¬ 
actions between the stars and the center of mass of the 
binary. As N rises the perturbations that the binary 
experiences become less noisy and the binary wanders 


less from the center. It therefore reduces the density 
near the center more and its hardening slows down. 
The wandering amplitude cannot be made arbitrarily 
small by raising A, as it could if the binary were a sin¬ 
gle particle, because the restoring force that keeps the 
binary at the center grows weaker as the central den¬ 
sity goes down. There is therefore an A value above 
which the hardening stops decreasing with A, about 
10 ® for the experiment considered here: panel (i) 
shows that the hardening rate with A = 2 x 10® is the 
same as with A = 10®. For smaller bhs the wander¬ 
ing continues growing to larger A values. The mass 
spectrum affects the wandering less than it affects the 
eccentricity, because the center of mass gets perturbed 
by all the stars near the center, not just by the stars 
having close encounters with the binary. The mean 
hardening rates found with equal-mass particles dif¬ 
fered only slightly (they were larger) from those found 
with a mass spectrum. 

This sample integration has thus shown that our 
program can integrate binaries accurately in large N- 
body systems without great expense, but also that the 
results have to be extrapolated with care before they 
can be applied to real galaxies. The value of A can 
easily be the largest source of error if it is small. For 
the eccentricity it is possible with the help of a mass 
spectrum to choose A large enough to suppress the 
noise—the value required is inversely proportional to 
the masses of the bhs. For the semimajor axis this 
is more difficult, because the mass spectrum does not 
help as much, but it is still possible if the bhs are large. 
The main uncertainty that the value of A introduces is 
in the hardening rate; it does not affect the response 
of the galaxy (except right at the center) if that is 
plotted versus the binary semimajor axis as it is in 
panel (c) of Figure ^ 

4 The effect of wandering on 
the hardening and ejection 
rates 

Our N-body program can be used both to test the pre¬ 
dictions from our 3-body scattering experiments and 
to study dynamical processes that could not be ex¬ 
plored with those experiments. In this section we test 
the predictions for the hardening and mass-ejection 
rates. There are several reasons for suspecting that 
these rates may differ in the N-body experiments. The 
most obvious is that the stellar distribution evolves 
along with the binary in the N-body experiments, 
whereas in the scattering experiments it is assumed 
to be fixed. When the low-angular-momentum parts 
of the distribution get depleted the hardening slows 
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Figure 2: Dependence of the binary evolution on the number of particles {N), for the experiment shown in Fig. 0. The 
left, middle, and right columns show the eccentricity and semimajor axis computed with N/W^ = 6.25, 25, and 100. 
The integrations in the top row used equal-mass particles (A = 0); those in the bottom two rows used a mass spectrum 
(A = 0.75). For each combination of A and N the integration was done three times with three different orbital planes 
for the initial orbits of the bhs. The dotted lines in panels (f) and (i) are from an integration with N = 2 x 10® and 
A = 0.75. The dashed lines in the bottom row are from panel (b) of Fig. 


down. But differences can occur even in the absence 
of loss-cone depletion, because in the N-body exper¬ 
iments the stars are influenced by forces from both 
the binary and the galaxy, not just from the binary 
as in the scattering experiments, and because in the 
N-body experiments the binary does not remain fixed 
in space. We argue here that the wandering of the 
binary increases the ejection rate by about a factor of 
two. 

The evidence for this is shown in Figure In 
panel (b) we plot the ejected mass as a function of 
the binary orbital velocity during integrations of bi¬ 
naries in a Jaffe model (with N = 10®, and with the 
binaries started as was the binary in Fig. We 
computed the ejected mass at time t from Mej{t) = 


M{rt:,0) — M{r^,t), with r* the radius at which 
p{rt,,t)/p{rt,,Q) = 0.95. The five lines in panel (b) 
are about the same within the uncertainties; we take 
the dashed-dotted line as the mean N-body result (in¬ 
tegrations with other N values gave similar results). 
We also plot the prediction from our 3-body scattering 
experiments, computed from 

/*00 J 

Mej(a)/Mi2 = C+ / — J(a), (8) 

where C, which we set to 0.6, is to allow for the mass 
that gets removed when the bhs first sink into the cen¬ 
ter, and J is the 3-body ejection rate (in what follows 
we speak of dM^-JdVhm loosely as though it were the 
ejection rate). In computing the 3-body prediction 
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Figure 3: Ejected mass during integrations of equal-mass binaries in a Jaffe model, (a) Stellar density at the start 
(top line) and end (bottom line) of the integration of Fig. |^. The dashed line is drawn at r = r*; the shaded region 
indicates the ejected mass, (b) Ejected mass versus binary orbital velocity from integrations like that of Fig. but 
with mi = m 2 = 0.04 (solid line), 0.02 (dotted), 0.01 (short dashed), 0.005 (long dashed), and 0.0025 (dashed dotted); 
the heavy solid line is the prediction from the scattering experiments of paper I. 


we assumed a Maxwellian velocity distribution with 
cr^ = 1/2, and, as we did in paper I, we counted stars 
as ejected if they get expelled from the binary with 
a velocity larger than Wej = V^cr. Although there is 
some uncertainty in the appropriate values for (7, cr, 
and Vej , it is too small to explain the factor-of-two dif¬ 
ference in the slopes of the N-body and 3-body ejec¬ 
tion curves near the end of the integration. Similar 
differences were found from integrations with galaxy 
models with other 7 values. 

A simple thought experiment suggests that the wan¬ 
dering of the binary can account for a large part of 
the difference. Suppose the binary moves so fast that 
incoming stars can interact with it only once before 
it moves out of their reach. For some stars this will 
not matter, but for others, which would experience 
multiple interactions before getting expelled if the bi¬ 
nary were not wandering, it will. The cross sections 
in Figure 4 of paper I show that the hardening rate 
is sensitive to these multiple scattering events. Addi¬ 
tional calculations from those cross sections show that 
the hardening rate for a hard, equal-mass binary gets 
reduced by a factor of 7.5 when stars are allowed to 


interact with the binary only once, whereas dMe^/dt 
gets reduced by at most a factor of two. This thought 
experiment is admittedly unrealistic, because in real¬ 
ity the binary will not be moving as fast as we have 
assumed, but the qualitative conclusion drawn from it 
should be correct: wandering will decrease the hard¬ 
ening rate more than dMej/dt, and hence will increase 
the ejection rate when it is expressed as dMej/dVhin 
or dMej/dln(l/a). 

To test this conclusion we compared the ejection 
rates from N-body integrations with the binaries free 
to wander with those from integrations with the cen¬ 
ter of mass of the binaries fixed at the origin (as was 
described in Sec. ||). We added a small triaxial per¬ 
turbation to the potential (the same for both sets of 
integrations) to allow more stars to reach the center, 
for otherwise in the integrations with the center of 
mass fixed the hardening would not have continued 
far enough for us to make a fair comparison. The 
ejection rate was close to the 3-body prediction when 
the center of mass was fixed, but not when it was free; 
the difference was again about a factor of two. 

If our explanation is correct, the hardening rate in 
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the N-body experiments should be lower than was pre¬ 
dicted from the scattering experiments. To test this 
we performed integrations of equal-mass binaries in a 
Plummer model, starting the binaries as did Mikkola 
and Valtonen (1992). The Plummer model is better 
than the Jaffe model for this test because its large, 
flat core simplifies the 3-body prediction. Mikkola 
and Valtonen said that their N-body hardening rate 
agreed with the 3-body prediction, but we found a 
lower hardening rate (we checked the dependence on 
all our integration parameters); their rate might have 
been too large by chance, since the results are noisy 
with only 10^ particles, the number that they used. 
To make a better comparison of the 3-body and N- 
body rates we performed a series of integrations like 
that of Mikkola and Valtonen but increasing the par¬ 
ticle number and decreasing the bh masses by a factor 
of two each time, i.e. using iV/10^ = 1, 2, 4, 8, 16, to¬ 
gether with nil = m 2 = 1/100, 1/200, 1/400, 1/800, 
1/1600. The importance of loss-cone depletion goes 
down as mi is lowered with Nmi held fixed, because 
the mass of the binary goes down but the wandering 
amplitude does not. We found the N-body harden¬ 
ing rate to be about 0.5-0.7 as large as the 3-body 
prediction. 

This change to the 3-body hardening rate is not 
important for estimates of merger times, given the 
much larger change that results from loss-cone deple¬ 
tion (we describe later how wandering changes the 
hardening rate when loss-cone depletion occurs). But 
the change to the ejection rate is important, because 
it means that before merging the binaries will eject 
about twice as much mass as we predicted in paper I. 
A question still to be answered is whether this is true 
for binaries with unequal masses. 

5 The factors affecting the ec¬ 
centricity evolution 

We now study the factors affecting the eccentricity 
of a binary during the early stages of the evolution. 
We first consider binaries with equal or nearly-equal 
masses, whose early evolution can be predicted by the 
usual dynamical-friction formula, and then binaries 
with a large mass ratio, whose eccentricity may be 
influenced by resonant interactions with the stars. 

5.1 Non-resonant friction 

Chandrasekhar’s dynamical-friction formula can be 
used in most circumstances to predict the loss rates for 
energy and angular momentum as a bh binary forms 
at the center of a galaxy. The eccentricity is influenced 
by two competing factors. The velocity dependence 


of the friction formula increases the frictional force at 
apocenter relative to that at pericenter if the bhs are 
moving fast, and therefore favors an increase in the ec¬ 
centricity (Fukushige et ah, 1992). The density gradi¬ 
ent in the galaxy has the opposite effect if the density 
is higher at pericenter than at apocenter (Polnarev 
and Rees, 1994). Anisotropy is a further complication. 
Casertano et al. (1987) integrated the orbits of single 
massive particles in isotropic and radially anisotropic 
galaxy models using both Chandrasekhar’s formula 
and N-body methods, and found that anisotropy can 
counteract the density gradient and allow orbits to 
reach the center without their eccentricity decreasing. 
The eccentricity did not increase, however, unless the 
anisotropy was strong enough to make the galaxy sus¬ 
ceptible to the radial-orbit instability. 

Since Chandrasekhar’s formula is applicable only 
to the early stages of the binary evolution, the pre¬ 
dictions that it makes for the eccentricity need to 
be checked by N-body experiments. Three previous 
studies have done this, all using a Plummer model 
to represent the galaxy. Makino et al. (1993) did 
experiments with 16,000 particles to study the for¬ 
mation and evolution of equal-mass binaries with 
bh masses ranging from 0.01 to 0.16. Their bina¬ 
ries sometimes formed with large eccentricities, but 
those experiments started with the bhs approach¬ 
ing the center on nearly radial orbits. Phinney and 
Villumsen (Phinney, 1994) did a similar experiment 
but started the bhs with a larger angular momentum 
(with V = 4 X 10"^, mi = 0.03, m 2 = 0.02) and found 
that the eccentricity remained small. Mikkola and 
Valtonen (1992) also started the bhs with a large an¬ 
gular momentum and found a small initial eccentric¬ 
ity, although their experiment started with the binary 
nearly formed, and the noise level was high enough 
(with N = 10^, mi = m 2 = 0.01) that the eccentric¬ 
ity bounced around in a random manner once the bi¬ 
nary began to harden. The results from these studies 
were thus consistent with the predictions from Chan¬ 
drasekhar’s formula, but were noisy and did not ex¬ 
plore the effect of a density cusp or a radial anisotropy. 

To improve upon this we performed experiments 
like those of Makino et al. (1993) but using a vari¬ 
ety of galaxy models and initial orbits, and using a 
large number of particles to suppress the noise. For 
the galaxy models we used those of Plummer, Hern- 
quist, and Jaffe, including for the Jaffe model both 
the isotropic model and a model with a constant ra¬ 
dial anisotropy of /3 = 1 — CT//cr^ = 0.5 (Cuddeford, 
1991). This [3 value is small enough to suppress the 
radial-orbit instability (Merritt and Aguilar, 1985), 
which would be absent from our calculations in any 
case because we include only spherical basis functions 
in the expansion of the potential. The bhs had masses 
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Figure 4: Eccentricity of an equal-mass binary when it first becomes hard at the center of a galaxy. The abscissa 
is the initial angular momentum of the bhs (at r = 0.5) in units of the angular momentum of a circular orbit of the 
same energy. The points and error bars show the mean eccentricity and the average deviation from the mean when the 
semimajor axis reaches a = 0.004, computed from 3-5 experiments with different orbital planes for the initial orbits. 
The isotropic (/3 = 0) Plummer and Hernquist models used N — 10® and A = 1; the Jaffe models used N = 2 x 10® 
and A = 0.5. 


mi = TO 2 = 0.01 and were started symmetrically 
about the origin at (r, v) and (—r, —v), with r = 0.5 
and with angular momenta equal to L/Lc times the 
angular momentum of a circular orbit of the same 
energy. The quantity that we measured is the eccen¬ 
tricity when the binary first becomes hard, i.e. when 
it reaches Uh = Gm^jAa^. 

The results in panel (a) of Figure ^ show that it 
is easier for a binary to form with a large eccen¬ 
tricity in the Plummer model than in the Hernquist 
model, in agreement with the prediction of Polnarev 
and Rees (1994) that a density cusp lowers the ec¬ 
centricity. The isotropic Jaffe model is similar to the 
Hernquist model this respect (the Jaffe-model results 
are noisier than the Hernquist-model results at small 
L/Lcj possibly because of the different velocity distri¬ 
butions in the two models). For the anisotropic Jaffe 
model the eccentricities are larger at small values of 
L/Lc but are still small if L/L(. > 0.6, in agreement 
with the prediction of Casertano et al. (1987) that a 
radial anisotropy can allow an eccentric orbit to re¬ 
main eccentric but cannot cause the eccentricity to 


grow. 

If the incoming bh orbits are distributed isotropi¬ 
cally, the probability of an orbit of a given energy hav¬ 
ing an angular momentum less than L is {L/Lc)'^. The 
probability of a binary forming with a large eccentric¬ 
ity is therefore small. The formation of a bh binary 
after the merger of two galaxies will be more compli¬ 
cated than in the simple experiments considered here: 
the BHS will have clusters of bound stars, and the sur¬ 
rounding galaxies will be greatly disturbed. Barnes 
(1992) simulated the merger of disk galaxies (with¬ 
out bhs) and found that the final collision of the two 
bulges was often remarkably head-on, suggesting that 
small values of L/Lc might not be that improbable. 
But, as Barnes admits, the bulges in his galaxy mod¬ 
els had unreasonably large core radii. If the bulges 
had been more concentrated their final collision might 
have had a larger L/Lc, for the same reason that the 
BHS in our calculations approach the center with a 
larger Lj when the galaxy has a density cusp. Thus, 
except for unusual cases in which a galaxy has a large, 
flat core or a strong radial anisotropy, we suspect that 
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most mergers of galaxies with central bhs will lead to 
binaries with small eccentricities. 

5.2 Resonant friction 

There is one circumstance in which dynamical fric¬ 
tion can cause the eccentricity of a bh binary to in¬ 
crease even in the presence of a density cusp. Rauch 
and Tremaine (1996) described this in their discus¬ 
sion of resonant relaxation, and called the frictional 
force “resonant friction” because it results from res¬ 
onant interactions between the stars and the smaller 
of the two BHS. Consider a bh of mass m 2 orbiting 
a larger bh of mass mi, with mi surrounded by a 
cluster of bound stars of total mass M^,, and assume 
that the stars have an anisotropic velocity distribution 
{df/dL ^ 0). If mi ^ max(m 2 ,M*), the stars move 
in a potential that is nearly Keplerian and have orbits 
that are nearly closed. The frictional force on m 2 then 
cannot be described by the usual dynamical-friction 
formula. Rauch and Tremaine show that resonant in¬ 
teractions between m 2 and the stars increase \dL/dt\ 
by a factor of order mi/max(m 2 ,M*) but have no 
effect on \dE/dt\, and thus cause the eccentricity of 
the binary to change much faster than the semimajor 
axis; the sign is such that the eccentricity increases if 
df/dL < 0 (i.e. if the distribution is biased in favor of 
radial orbits) and decreases if df/dL > 0. 

To determine the strength of the anisotropy re¬ 
quired for resonant friction to be effective, we made 
five galaxy models having the same mass distribution 
but different anisotropies, each model containing a 
central bh of mass mi = 0.1 (see Fig. ||). We chose a 
large value for mi so that we could choose m 2 ^ mi 
and still have m 2 much more massive than the stars; 
one can think of our galaxy models in these experi¬ 
ments as representing the inner parts of real galaxies. 
Model D was made by growing mi slowly at the center 
of a 7 model with 7 = 0, which should lead to a final 
model with a cusp p ~ and a central anisotropy 
(3 ~ —0.2 (Quinlan et ah, 1995); the central /3 value in 
the hgure is slightly lower than the theoretical predic¬ 
tion because we removed some of the stars on radial 
orbits close to the bh. The other models were made 
by varying the cusp slope 7 of the initial model and 
the growth time of the bh (short growth times lead 
to radial anisotropies, long growth times to tangen¬ 
tial anisotropies), and sometimes by removing stars 
on orbits with low or high angular momenta after the 
BH growth was finished to get the desired mass and 
anisotropy profiles. The models were integrated for 
a while to check that they had reached equilibrium. 
Resonant friction should be important when a small 
mass m 2 sinks inside r < 0.01, where M{r) <C mi and 
where the models have different anisotropies. 


The orbits of small bhs were integrated in the five 
models starting at r = 0.2 and continuing until the bi¬ 
nary semimajor axis reached a < 0.001. To reduce the 
cpu time, mi was fixed at the origin, which is what 
Rauch and Tremaine assumed. This assumption is 
not valid near the end of the integrations, when the 
binary is hard and the motion of mi would be impor¬ 
tant, but it does not affect the eccentricity much when 
m 2 first enters the region r < 0.01; we checked that 
by repeating a few of the integrations without fixing 
mi. Figure ^ shows the eccentricity evolution for m 2 
values 0.001, 0.003, and 0.01. The eccentricity and 
semimajor axis in this figure are generalized orbital 
elements that agree with the Kepler elements when 
the BHS are close and give useful information when 
they are not close (the semimajor axis then gives an 
estimate of the distance between the bhs). Since res¬ 
onant friction should increase the eccentricity if /3 > 0 
and decrease it if /3 < 0, the initial eccentricities were 
chosen to be small for the radial models (/3 > 0) and 
large for the tangential models (/? < 0). 

Consider first the results in panel (a) from the 
integrations with m 2 = 0.001. The influence of 
anisotropy on the eccentricity is clear: in models A 
and B, with radial anisotropies, the eccentricity grows 
larger, whereas in models D and E, with tangential 
anisotropies, it grows smaller. In models A and E 
the eccentricity changes much faster than the semi¬ 
major axis when 1/a ~ 100, as Rauch and Tremaine 
predicted. The evidence from the integrations with 
larger m 2 values is not as clear. The eccentricity still 
grows in model A when m 2 = 0.003, suggesting that 
resonant friction still operates, but the growth is not 
as fast as with m 2 = 0.001. When m 2 = 0.01 the 
eccentricity goes down in all the models, and the rate 
of change is low enough to be explained by ordinary, 
non-resonant friction. This m 2 value perturbs the res¬ 
onances responsible for resonant friction. 

Thus resonant friction can cause a rapid growth in 
the eccentricity of a bh binary, but only if two con¬ 
ditions are satished: first, m 2 must be at least 30 
times smaller than mi (this was implicit in the for¬ 
mulas of Rauch and Tremaine); second, and more 
troublesome, the cluster of stars bound to mi must 
have a strong radial anisotropy, say /3 > 0.5. There 
is no reason to expect an anisotropy like that around 
a massive bh. The only existing models for the for¬ 
mation of BHS in galaxies that make definite predic¬ 
tions for the anisotropy—the adiabatic-growth model 
(Young, 1980; Goodman and Binney, 1984; Quinlan 
et ah, 1995), and the binary-merger model studied 
here—predict a tangential anisotropy. Any model in 
which the bh grows by consuming stars on radial or¬ 
bits is likely to favor a tangential anisotropy. Some 
formation models may allow a radial anisotropy (e.g.. 
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Figure 5: Anisotropy /3 = 1 — a^ja^ and enclosed stellar mass M for the inner parts of five models nsed to study 
resonant friction. Each model has N = 50K particles with a mass exponent A = 0.75, and a central bh of mass 
mi = 0.1 with a BH-star softening length of 0.00025. 


models that involve a sudden collapse at the center, 
with the BH and the stars forming simultaneously), 
but they have not been studied and are unlikely to pre¬ 
dict an anisotropy as strong as /3 > 0.5. There are two 
galaxies for which we have estimates of the anisotropy 
around the central BH. If the BH mass estimate of 
Harms et al. (1994) for M87 is correct, then the best 
fit to the stellar velocity distribution has a tangential 
anisotropy, although the error bars are large enough 
to allow a slight radial anisotropy (Dressier and Rich- 
stone, 1990; Merritt and Oh, 1996). For M32 also the 
best fit to the velocity distribution around the central 
BH has a tangential anisotropy (van der Marel et al., 
1997). 

Although we do not believe that resonant friction 
will be important for the eccentricity evolution of bh 
binaries, the possibility of this is nevertheless inter¬ 
esting: it shows the danger of relying on the simple 
dynamical-friction formula. Resonant friction may be 
more important for binaries in galaxies with disks, 
where it will affect the inclination as well as the ec¬ 
centricity. 


6 The response of the galaxy 

The response of a galaxy to the hardening of a BH 
binary may explain why large elliptical galaxies have 
weaker density cusps than smaller galaxies. Figures 
|] and H showed some results for the core expansion 
and mass ejection caused by a binary at the center 
of a Jaffe model. Here we examine the change in the 
density proHle for a range of binary masses, and also 
the kinematical signature that the binary leaves in the 
velocity distribution. We again use a Jaffe model to 
represent the galaxy; the results will apply with some 
changes to models with other gamma values. 

We integrated binaries with hve masses (mi = 
m 2 = 0.04, 0.02, 0.01, 0.005, and 0.0025), starting 
them as we did for the integration in Figure |^. The 
galaxy model had 10^ particles with a mass exponent 
A = 0.75. To study the change in the galaxy we kept 
frequent records of the radii containing fixed fractions 
of the total mass, like the radii in panel (c) of Fig¬ 
ure 1^ but more of them, and also of the sums of mv^ 
and mU( for the stars inside those radii. By differ¬ 
entiating the mass fractions and velocity sums with 
respect to radius, after first smoothing their time de¬ 
pendence, we were able to recover the density and the 
radial and tangential dispersions at any time during 
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Figure 6: Eccentricity versus semimajor axis for bhs of mass m 2 = 0.001 (a), 0.003 (b), and 0.01 (c) orbiting in the 
five models of Figure ^ (the different line types correspond to those used in Fig The two lines for each model and 
m 2 value are from integrations in which m 2 was started in a clockwise and counter-clockwise sense. Panel (d) shows 
five of the lines from panel (a) before they were smoothed. The orbital elements are dehned (for this figure only) by 
1/a = —2E/mi and e = yl —with E and L the energy and angular momentum of m 2 and with the potential 
of the galaxy set to zero at r = 0.2 for the calculation of E. 


the integrations. We use this data in Figure to plot 
the densities and dispersions when the five binaries 
reached Vbin = 4, the highest velocity reached by the 
binary with toi = 0.04 (we stopped the integration 
then because it was becoming slow). The binaries 
with smaller masses reached higher velocities, but to 
study the dependence on the binary mass it is best to 
compare the galaxies at the same Vbin- 

The density profiles in panel (a) show that a binary 
can eliminate a density cusp. Although the profile 
for the smallest binary has a mild cusp, the others 
have densities that are flat or, in some cases, decreas¬ 
ing towards the center; they would probably all de¬ 
crease if we could plot them closer to the center (this 
was predicted by Begelman et al. 1980). The mass 
ejected from the central region is Mej/Mi 2 = 1.2-1.5, 
as we reported in Figure H, with no systematic de¬ 
pendence of Mej/Mi 2 on the binary mass. Because 
the Jaffe model has a mass that, near the center, in¬ 


creases linearly with radius, the core radius at a fixed 
Vbin increases linearly with the binary mass; for other 
gamma models the dependence would be nonlinear. 
The Mej/Mi 2 values computed here are larger than 
we predicted in paper I, partly because some mass 
gets displaced from the center in the N-body exper¬ 
iments before the binaries form and partly because, 
as we discussed in Section |[ the ejection rates once 
the binaries form in the N-body experiments (and, we 
believe, in real galaxies) are higher than we predicted 
in paper 1. As the binaries continue hardening, more 
mass will be ejected and the core radii will grow. For 
example, Figure ^ shows the density profile for the bi¬ 
nary with mi = 0.01 at Vbin = 6.8: the central density 
is less than half of its value at Vbin = 4, and the core 
radius is nearly twice as large. The final profile will 
depend on how far the binary has to harden before 
gravitational radiation becomes important. Massive 
BH binaries in typical galaxy cores have to shrink by 
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Figure 7: Response of a Jaffe model to the hardening of a bh binary. The dashed lines show the initial (a) density, 
(b) 3D velocity dispersion, and (clanisotropy. The solid lines show these quantities when the binary reaches Vbin = 4, 
with the binary started as in Fig. hi and with the five lines for bh masses mi = m 2 = 0.04 (rightmost line), 0.02, 0.01, 
0.005, and 0.0025 (leftmost line). 


a factor of 10-50 beyond the point at which they be¬ 
come hard, with the factor varying with the bh masses 
as for an equal-mass binary. The binaries in 

Figure 0 shrank about a factor of 4 beyond the point 
at which they became hard. In a typical galaxy they 
would have to shrink another factor of 3-12, causing 
I4in to rise by a factor of 1.7-3.5 and, using Figure || 
as a guide, adding 0.4-0.8 to the value of Mej/Mi 2 . 

The kinematical response of the galaxy is shown 
in panels (b) and (c). The velocity dispersions rise 
towards the center, approximately (but not exactly) 
as ~ l/r. Mass estimates derived from these dis¬ 
persions will underestimate the central masses if the 
anisotropies are not taken into account. The mod¬ 
els all develop strong tangential anisotropies, because 
the stars on radial orbits are the ones that get ejected. 
Four of the five lines in panel (c) show /3 ~ — 1 near 
the center, although we see no reason why that value 
should be preferred. The central anisotropy is difficult 
to measure (the leftmost parts of our lines are affected 
by noise); we expect that (3 would fall below —1 if we 


could measure it closer to the binary. This is a much 
stronger anisotropy than is predicted by models like 
the adiabatic-growth model for the formation of mas¬ 
sive black holes. The adiabatic growth of a bh in an 
initially isotropic galaxy leads to a final anisotropy of 
at most (3 = —0.3 (Quinlan et ah, 1995) 

Our mass-ejection results support the suggestion of 
Ebisuzaki et al. (1991) and Faber et al. (1996) that bh 
binaries can explain the weak density cusps in large el¬ 
liptical galaxies. If these galaxies contain massive bhs, 
and if they have formed through mergers of smaller 
galaxies also containing massive bhs, the binaries that 
result will clear out the central region. A real galaxy 
merger will of course be more complicated than the 
idealized experiments considered here. In most cases 
a binary will result when a small satellite galaxy sinks 
into a large galaxy containing a massive bh. The scat¬ 
tering experiments from paper I suggested that more 
mass gets ejected if a massive bh merges with a num¬ 
ber of small BHS than if it merges with one large bh 
(with the same total mass), although each small bh 
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may counteract the mass ejection to some extent by 
bringing fresh stars and gas to the center. The general 
result from our experiments should remain true: the 
density cusp in a merged galaxy will be weaker if the 
merger leads to a bh binary. And with this result we 
can make a testable prediction: if this is the explana¬ 
tion for the weak density cusps in large galaxies, the 
central velocity distributions in those galaxies should 
have strong tangential anisotropies. As we mentioned 
in Section there is evidence that the galaxy M87 
passes this test. 

7 The binary merger time in 
real galaxies 

The answer to many questions involving massive bh 
binaries depends on the binary merger time. A big 
uncertainty in computing this has been loss-cone de¬ 
pletion. Although the loss cone can be replenished if 
the potential is nonspherical, and the binary merger 
can be helped by gas accretion onto the BHs, for 
many galaxies these perturbations will not be enough. 
We therefore consider whether a binary can merge 
through purely stellar dynamical processes in a spher¬ 
ical galaxy. We believe that the answer is often yes, 
thanks to the wandering of the binary from the galaxy 
center. 

Previous discussions of the loss-cone problem have 
mostly assumed that the central object (a single bh 
or a binary) remains fixed at r = 0. The rate at 
which stars accrete onto or interact with the object 
once the loss cone becomes depleted is then given by 
an expression of the form (e.g., Shapiro, 1985) 

dM ^ p{r„it)r^rit 
dt tr(^crit) lB(rcrit/^D) 

where td is the destruction or interaction radius, 
which we take to be the binary semimajor axis, and 
Tcrit is the radius at which the loss-cone angle 9\c 
equals the rms deflection suffered by stars through 
two-body relaxation in one dynamical time. For 
BH binaries in typical galaxies, rcrit is larger than 
GM 12 /, and the relaxation time at rent is much 
longer than the age of the galaxy. The hardening rate 
implied by equation (||) is then too low to allow most 
binaries to merge. 

Valtonen (1996) used the galaxy M87 to illustrate 
how difficult it is for a central binary to merge when 
the relaxation time is long. He considered an equal- 
mass binary with M 12 = 3 x 10® Mq and a separation 
a = 0.6 pc, at which the merger time from gravita¬ 
tional radiation is 10 ^® yr, and at which the total mass 
in stars that can approach within a of the center— 
assuming that the loss cone is not depleted—is only 


1/200 of the mass of the binary. For the binary to 
shrink by a factor of e, the loss cone would have to 
be refilled about 200 times. The problem is worse for 
smaller bhs because they have to shrink more before 
radiation becomes important. For M 12 = 3 x 10® Mq 
the factor of 200 in Valtonen’s example gets replaced 
by 7800. Valtonen concluded that, for nearly the 
whole range of parameters of interest, bh binaries in 
galaxies are incapable of merging. 

Wandering changes this conclusion. We refer here 
not to the small correction to the hardening rate that 
we discussed in Section |[ but to a larger correction 
that occurs when the loss cone becomes depleted. A 
simple equipartition estimate predicts that a particle 
of mass Mi 2 (a single BH or a binary) in a homoge¬ 
neous core of smaller particles of mass m will have a 
maximum wandering amplitude 
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Young (1977) argued that this amplitude is sometimes 
large enough to affect the growth of bhs in the cen¬ 
ters of galaxies, but most subsequent work has either 
ignored wandering or has dismissed it as unimpor¬ 
tant. For BH binaries we believe that wandering is 
even more important than Young suggested, because 
the prediction (0) underestimates the wandering am¬ 
plitude of a BH binary. The prediction is applicable 
to a massive particle in a homogeneous core, but a 
BH binary ejects stars from the center of the core, re¬ 
ducing the restoring force that keeps it at the center, 
and the superelastic encounters between the binary 
and the stars invalidate the equipartition assumption 
on which the prediction is based. Consider the inte¬ 
gration that we used earlier as an example. The core 
radius in Figure ^ is Tc ~ 0.04 (derived by fitting a 
King-model core to the final density), and the mass 
ratio is m/Mi 2 < 10 “®/ 0.02 (this would be an equality 
if the stars had equal masses), so the prediction ( 10 ) 
gives Tw ^ 0.0009. But the maximum wandering am¬ 
plitude shown in Figure ^ is more than five times as 
large. And, as we mentioned in Section the ampli¬ 
tude remained the same when we raised N from 10® 
to 2 X 10®, making the discrepancy with the prediction 
even larger. 

We suspect that in many galaxies it is the wander¬ 
ing of the binary, and not diffusion of stars into the 
loss cone, that determines the hardening rate once the 
loss cone gets depleted. In our Jaffe-model integra¬ 
tion, for example, the hardening rate did not increase 
when we added non-spherical terms to the potential 
by changing (max from 0 to 2 , even though this change 
made the potential noiser, somewhat like the way that 
two-body relaxation does, and allowed more stars to 
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reach the center. The change did increase the harden¬ 
ing rate when we applied a constraint force to fix the 
binary at the center of the galaxy. 

Wandering will not be enough to allow all bh bina¬ 
ries to merge. It allows a binary to continue hardening 
beyond the point at which it would stop without wan¬ 
dering, but the hardening continues at a reduced rate. 
As a rough guide, we can say from our integrations 
that once a binary ejects the stars from its immediate 
vicinity and creates a core, it hardens at a rate that is 
about 10-50 times slower than the rate we predicted 
in paper I ignoring loss-cone depletion. Whether this 
will be fast enough to allow a merger depends on the 
particular galaxy under study. Gas and nonspherical 
perturbations will help if they are present. 

These conclusions have been derived for equal-mass 
binaries. The most likely formation scenario for a 
massive bh binary is for a small satellite galaxy con¬ 
taining a central bh to sink into a larger galaxy with a 
larger bh. The bottleneck in this case may be the time 
required for the small bh to reach the center. N-body 
experiments by Balcells and Quinn (1989, 1990) show 
that small, dense satellites can easily sink into the cen¬ 
ter of large elliptical galaxies, but Weinberg’s (1997) 
perturbation calculations suggest that this is not al¬ 
ways the case. If the satellite galaxy surrounding the 
small BH gets disrupted, the bh will take longer to 
reach the center. Neither Weinberg nor Balcells and 
Quinn included bhs in their models; some preliminary 
experiments with bhs have been done by Governato et 
al. (1994). More work on this is needed, since the pre¬ 
dicted merger rate for the gravitational-wave detector 
LISA is much larger if small bhs can participate in 
the mergers (Haehnelt, 1994). 
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A The hybrid N-body program 

We describe here some technical details of our hybrid 
N-body program. We assume that the star cluster 
contains two bhs, but the program can work with any 
number of bhs. Our notation follows that of Aarseth 
(1994) and Hernquist and Ostriker (1992). 


A.l The integrator 

The program uses different integrators for the bhs and 
the stars. The stars are integrated with the NBODYl 
integrator. The stepsizes for the stars vary continu¬ 
ously and are updated after each step by a function in¬ 
volving the force and the first three time derivatives of 
the force, multiplied by the square root of an accuracy 
parameter ry. We used rj = 0.005-0.01, smaller than 
the value that Aarseth recommends {rj = 0.02-0.03), 
so that we could integrate close encounters accurately 
with a small softening length. The bhs are integrated 
with the NBODY2 integrator, which uses the Ahmad- 
Gohen neighbor scheme. The bhs have two stepsizes, 
6ti and Jfo, determined by two accuracy parameters, 
rji = rj and rjr = 2rji. The bhs are moved and the irreg¬ 
ular force from nearby stars (the neighbors) is updated 
on the smaller stepsize 5ti] the regular force from the 
remaining stars is updated on the larger stepsize 5tr 
and is predicted between those times from its known 
time derivatives. The number of neighbors for each 
BH is set to approximately (A^/4)^/^, the scaling rec¬ 
ommended by Makino and Hut (1988); the program 
forces each bh to be in the neighbor list of the other 
BH at all times. When N is large, 5ti is usually much 
smaller than Str- The neighbor scheme improved the 
efficiency of our integrations by at least a factor of 
ten. 

The program checks periodically to see if the bhs 
are close enough for their motion to be regularized. 
If they are, the bhs are replaced by a center-of-mass 
particle and a set of KS coordinates describing their 
relative motion. Our criterion for “close enough” is 
that each bh is the closest neighbor of the other bh, 
that the force on each bh from the other bh is five 
times larger than the force from the stars, and that 
the neighbor list for the center-of-mass particle is large 
enough to contain the perturbers for the KS coor¬ 
dinates, which means the stars that exert a dimen¬ 
sionless perturbation exceeding a value 7min- The KS 
equations of motion are integrated with the Hermite 
integrator from the NBODY6 program; the center-of- 
mass motion is integrated with the NBODY2 integra¬ 
tor as is the motion of single bhs, except that the 
forces between the center of mass and nearby stars 
(those within a distance Xa) are computed by resolv¬ 
ing the binary into its two bhs. We used A = 1000, 
probably much larger than was necessary (Aarseth 
recommends A = 70). The KS stepsize is chosen so 
that the number of steps per unperturbed period is 
27r/?7„, with rju the KS accuracy parameter (Aarseth 
recommends rju ^ 0.1; we used 0.07). If the binary 
has no perturbers, the Keplerian motion is integrated 
exactly. Even with perturbers, regularization allows 
close binaries to be integrated efficiently regardless of 
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their eccentricity, with no softening needed for the bh- 
BH interaction. 

The large difference between the masses of the bhs 
and the stars in our experiments caused some diffi¬ 
culties with the standard KS integration procedures. 
When we used the value for 7min that Aarseth rec¬ 
ommends (10“®) the stepsize for the center-of-mass 
particle was sometimes much smaller than the step- 
size for the KS coordinates, because a small star near 
the massive binary does not get included in the list 
of KS perturbers if its dimensionless perturbation is 
less than 7minj but it does get included in the center- 
of-mass neighbor list and it affects the center-of-mass 
stepsize. Large errors can result if there are many such 
stars. To improve the accuracy we reduced 7min to 
10“^, and we introduced a parameter T = 50 so that 
stars within a distance To of the binary are included 
in the KS perturber list whatever the size of their 
perturbation. Although these precautions slowed our 
integrations somewhat, the integration of close bina¬ 
ries still went several times faster when regularization 
was used. 

Because the program regularizes only the motion of 
the BH binary, and not the motion of stars approach¬ 
ing the binary, the interactions between the BHs and 
the stars have to be softened; otherwise some stars 
would be given tiny stepsizes when they encounter the 
BHS. We chose the softening length so that the soft¬ 
ening velocity (eq. was at least 2-3 times as large 
as the velocity of the binary at the end of the in¬ 
tegration; three-body scattering experiments showed 
that a softening length of that size does not affect 
the average energy transfer. Even with softening we 
sometimes had trouble with stars that got captured 
into tightly bound orbits around the bhs. We there¬ 
fore merged stars with the bhs if they were captured 
into orbits with semimajor axes less than twice the 
softening length. The number of stars merged in this 
way was small: the growth of the bh masses was at 
most 0.2%. 

A.2 The expansion of the potential 

The coefficients Anim for the potential expansion are 
updated at fixed time intervals At, after the coor¬ 
dinates of the stars are first predicted to a common 
time. We used the basis functions of Hernquist and 
Ostriker (the lowest-order member of the set is the 
Hernquist potential), but adjusted the scale length d 
for each initial model to customize the fit (Hernquist 
and Ostriker assume that d = 1). We found it nec¬ 
essary to soften the basis functions by replacing the 
radius r of a particle by (r^ -I- with h ~ 10“"*, 

for otherwise our use of both a mass spectrum and an 
integrator with individual stepsizes sometimes caused 


stars to get trapped in orbits with tiny stepsizes at 
the center of the galaxy. In most of our integrations 
we used a spherical basis (Zmax = 0), with the number 
of functions in the range 10 < rimax A 20, depending 
on the initial model (the Jaffe model required largest 
?T-max)- With a spherical basis the expansion assumes 
that the center of mass of the galaxy is fixed at the ori¬ 
gin, an assumption that has been shown to artificially 
accelerate the orbital decay of satellite galaxies in the 
sinking-satellite problem (White, 1983), although only 
when the satellite is far from the center of the primary 
galaxy (Quinn and Goodman, 1996). The assumption 
should not affect the evolution rate in our BH-binary 
experiments, because the bhs start close to the cen¬ 
ter and often symmetrically about the center, making 
the dipole distortion of the galaxy small. When we 
used a spherical basis the cpu time spent computing 
the expansion functions was at most half of the total 
cpu time; most of the time was spent computing the 
forces between the bhs and the stars. 

The NBODYl integrator needs the forces and the 
first three time derivatives of the forces acting on the 
stars at the start of an integration. This initialization 
is cumbersome when the expansion method is used. 
In principle the derivatives could be computed ana¬ 
lytically from the expansion functions, but in practice 
we used the known derivatives for the initial galaxy 
model, or used numerical interpolation to approxi¬ 
mate the derivatives. This would have been simpler 
if we had used the Hermite integrator for the stars, 
as that requires only the first time derivative (Makino 
1991). 

A.3 Tests and possible improvements 

We ran a number of tests to check the integrators and 
the potential expansion. We first fixed the expansion 
coefficients and checked that the conservation of en¬ 
ergy improved in the manner expected for Aarseth’s 
fourth-order integrators when we reduced the accu¬ 
racy parameters f] and rjy,. We also checked that the 
evolution of a hard binary was independent of which 
integrator we used for the bhs (NBODYl, NBODY2, 
or NBODY6), provided that the accuracy parameters 
were small (the required cpu time, of course, depends 
on the integrator). To test the potential expansion we 
checked that the program could maintain equilibrium 
7 models without bhs for radii r>10“^ (the minimum 
radius depends on 7 and Umax)- We also repeated 
some of our BH-binary experiments using the basis 
functions of Glutton-Brock (1973) instead of those of 
Hernquist and Ostriker, to check that the two gave 
the same result when Umax was large. 

One test did reveal an important limitation of our 
program. We ran some cold-collapse experiments 
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(starting a galaxy with the stars at rest, with no bhs) 
to check that the SCFBDY and SCF programs gave 
the same result for the virial ratio versus time. They 
did, but only if we used a stepsize At with SCFBDY 
considerably smaller than the stepsize required with 
SCF. The reason is that with SCFBDY the updat¬ 
ing of the expansion coefficients is not time reversible. 
The energy error after a fixed time therefore varies lin¬ 
early with At with SCFBDY, and not quadratically 
as it does with SCF. This did not matter for our work, 
because the change in the potential was slow and we 
were able to choose At small enough to make the er¬ 
rors small, but it would matter for a problem with a 
rapidly changing potential. 

Some improvements to the program may be pos¬ 
sible. The particles could be given block stepsizes 
to make the program run faster on vector comput¬ 
ers, although the vectorization will always be hindered 
by the complicated integration algorithm. The need 
to soften the BH-star interactions could be removed 
by more-sophisticated regularization, although this 
would be difficult when many stars interact with the 
BBS simultaneously. And the limitations of the poten¬ 
tial expansion—such as the linear dependence of the 
error on At, and the difficulty of resolving small stel¬ 
lar clumps away from the center—could be reduced 
by updating the contribution to the expansion coef¬ 
ficients from the inner part of the galaxy more fre¬ 
quently than the contribution from the outer part, or 
by using the direct-summation method for the inner 
part and an expansion method only for the outer part, 
although either of these changes would increase the 
program complexity, and using the direct-summation 
method would increase the cpu time. We did not try 
any of these ideas. 
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